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Abstract 

A novel scheme is presented for the direct numerical simulation of two-way coupling between a particle and 
an incompressible fluid flow maintained in a high-Reynolds-number turbulent regime. The main idea consists 
in combining a Fourier pseudo-spectral method for the fluid with an immersed-boundary technique to impose 
the no-slip boundary condition on the surface of the particle. This scheme is shown to converge as the power 
3/2 of the spatial resolution. This behavior is explained by the L2 convergence of the Fourier representation 
of a velocity field displaying discontinuities of its derivative. Benchmarking of the code is performed by 
measuring the drag and lift coefficients and the torque-free rotation rate of a spherical particle in various 
configurations of the carrier flow. Such studies show a good agreement with experimental and numerical 
measurements from other groups, and validate the code. A study of the turbulent wake downstream the 
sphere is also reported. The mean velocity deficit is shown to behave as the inverse of the distance from the 
particle, as predicted from classical similarity analysis. This law is reinterpreted in terms of the principle of 
"permanence of large eddies" that relates infrared asymptotic self-similarity to the law of decay of energy 
in homogeneous turbulence. 
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1. Introduction 



The dynamics of suspended particles whose size is larger than the smallest active scales of the carrier 
turbulent flow is important for the understanding of many natural and industrial applications. Let us 
mention two concrete problems borrowed from the study of planet formation and of rain processes. In the 
first situation, the role of turbulence is still unsettled in the d ynamical mechanisms t hat t ake place during 
the origination of planets in a proto-planetary disk (see, e.g.. Ide Pater and Lissauerl . 120011 ). A key issue is 
to understand whether or not the disk turbulence can generate an increase of the rate at which a large body 
can accrete dust. In other terms, it requires identifying the turbulent fluctuation mechanisms that lead 
to a snowball effect where the larger bodies have a growth speed all the more fast as their size increases. 
The same kind of question arises in the precipitation process when a large raindrop falls under the effect 
of gravity and collides with much smaller droplets that are transported by the turbulent airflow in the 
cloud (see IShaw , 2003 ) . For this problem, another question is to know the effi ciency of wet deposition (o r 
scavenging) of particle pollutants that are present in the atmosphere (see, e.g.. ISeinfeld and PandisL 1998 ). 
These various phenomena are still poorly understood and are subject to empirical ad-hoc modeling. Such 
fundamental open questions, which are taken from natural sciences but also arise in a number of industrial 
settings, reflect an inadequate current level of modeling for large particles in turbulent flows. 

Current approaches for the dynamics of finite-size particles are commonly assuming that the particles 
are infinitely small (point particles) and that the Reynolds number associated to their interaction with the 
flow vanishes, so that their dynamics can be obtained by matching the solution to the Stokes equation 
in their vicinity to an unperturbed fluid velocity. Recent exp erimental work showed that the dynamics 
of finite-size par t icles is not well described by such models (see Qureshi et al. . 2007 ; Xu and Bodenschatzi 
120081: IVolk etahl . l2010h . Significant size effects have been for instance measured in the acceleration statistics 
of such particles. Despite the fact that the flatness of the acceleration distribution depends only weakly 
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on the diameter of the particle, its variance strongly decreases with diameter and temporal correlations of 
the acceleration vary with the particle size in a manner that cannot be predicted by standard point-particle 
models. Improved models accounting for the so-called "Faxen corrections" (spatial averages of fluid velocities 
and accelerations over the volume of the particle) are in much better agree ment with such qualitativ e trends 
but hardly reproduce quantitatively the mentioned experimental data (see ICalzavarini et al. . 120091) . 

One of the most challenging requirement for improving such models is to get a precise understanding on 
how the carrier flow is modified in the vicinity of the particle in situations where such effects occurs on scales 
which are comparable to the active scales of the surrounding turbulence. For instance, it is clear that much 
work is still needed in order to get a better handling on how the kinetic turbulent energy, the dissipation rate, 
as well as other standard statistical properties of turbulence, varies in the turbulent boundary layer enclosing 
the particle. Direct numerical simulations appear to be valuable candidates as they give a unique global and 
instantaneous access to such quantities. However, the numerical problem is challenging for several reasons: 
(i) the numerical approach has to resolve all active scales of a fully developed turbulent flow, which span 
several decades; (ii) at small scales the boundary layer around the particle is of great importance and has 
to be well resolved; (Hi) if one considers a moving particle the algorithm has to follow the displacement of 
the particle accurately and efficiently; (iv) for a sufficient separation of the different scales involved (forcing 
scales, inertial range of scales, dissipative scales) the algorithm has to parallelize efficiently on massively 
parallel computers. 

Two different strategies have been recently proposed in order to solve the incompressible Navier-Stokes 
equations together with no-slip boundary conditions at the surface of a spherical particle. Both make use of 
a homogeneous grid covering of the complete physical do main. The first, "Physa l is", w hich was developed by 
Prosperetti and Qguz ( 2001 ) and extended and used bv lNaso and Prosperetti (2010), makes use of known 



analytic solutions to the Stokes equation for the boundary layer in the vicinity of the particle, in order 
to describe the velocity at the nearest surface points of the grid. Iteratively, the analytic solution and 
the solution on the homogeneous grid are matched. An advantage of this strategy is that the no-slip 
condition is imposed in an exact manner. However, improvement are required in order to achieve very-large 
R eynolds number carrier flows (the value Re\ — 20 of the Taylor-microscale Reynolds number was chosen 
in Naso and Prosperetti . 2010). The other recently developed approach consists in using and matching 



two different grids. A spherical one encloses the particle surface and is designed to resolve the boundary 
lay er and a second homogen eous grid covers the entire domain. This "Overset Grid" technique proposed 
by iBurton and Eatonl (|2005t) uses a third-order interpolation to transfer values from one grid to the other 
in order to couple the two solutions. This method was shown to reproduce very accurately experimental 
measurements of the drag coefficient at a particle Reynolds number of Re p = 20. In turbulent settings, this 
method was used with a Reynolds number Re\ = 32. 

Both of these methods reside on finite-difference schemes and are at the moment limited to rather small 
carrier-flow Reynolds numbers. We propose here a new method that allows one to reach much higher 
Reynolds numbers. The idea is to combine a standard pseudo-Fourier-spectral with a penalty method. The 
former is well adapted to incompressible homogeneous and isotropic turbulence and has been extensively 
used in such settings during the past. Its advantages comprise its degree of accuracy and its performance 
on massive parallel supercomputers such as the BlueGene/P system or large Intel/AMD clusters. In combi- 
nation with an immersed boundary and penalization strategy, high Reynolds numbers turbulent flow with 
an extended inert i al ran ge of scales can be reached. The proposed method is a modification of that used by 



Homann and Bed (2010) for investigating the dynamics of neutrally buoyant particles suspended in turbu- 



lent flows. The former version of the method was well-adapted to particles displaying a small slip velocity 
with respect to the fluid but happened to show some limitations when the slip is too large. We present 
here a modification of this method that unravels such shortcomings. This paper is organized as follows: 
in Section [2] we present in details the proposed method, together with some considerations concerning its 
convergence. Section [3] presents measurements of drag and lift coefficients in various settings and compare 
them to other numerical and experimental works. We consider the cases of a stationary sphere in a flow 
at a constant speed, of a freely rotating sphere and a translating sphere in a flow maintained at rest. Sec- 
tion ideals with the characterization of the turbulence decay in the particle wake. Section [5] encompasses 
concluding remarks and open questions. 
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2. The numerical method 



The dynamics of an incompressible flow in which is embedded a particle is given by the Navier-Stokes 
equations 

d t u = C{u) = -u-Vu-Vp + vV 2 u + f e , V • u = 0, for x E il\n p (1) 
together with the no-slip boundary condition on the surface of the particle 

u(x,t) = v p (x,t) for x 6 dQ p (t). (2) 

f e denotes here external forces that are exerted on the flow, fi is the full physical domain containing both 
the fluid and the particle, cT2 p (t) is the boundary of the physical location f2 p (t) C 51 of the particle at time 
t, and v p denotes the velocity of the particle surface. Note that in the case of a solid spherical particle with 
diameter d, one has f2 p = {x : \x — X p (t)\ < d/2} and v p (x,t) = V p (t) + [x — X p (t)] x fl p (t), where V p 
and flp are the translational and angular velocities of the particle, respectively. The particle motion is then 
obtained by solving Newton equations where the force exerted by the fluid is obtained by integrating on the 
surface of the particle the contribution from pressure and viscous stress to the fluid stress tensor. 

The numerical method presented here consists in integrating the Navier-Stokes equation ([!] ) by using 



a stro ngly-stable, low-storage Runge-Kutta scheme of third order that was introduced by and Oshcr 



(1988). It reads 



u n+i _ l u n + ! u » + £ At £(«"), where u" = ~u n + -«' + -At £(«') and u' = u n + At £(u°). (3) 
333 444 

At denotes here the time step and u n and u n+1 are the velocity fields at time steps n and n + 1 and the 
operator £ denotes the right-hand side (rhs) of (fT]). The later is integrated on a homogeneous rectangular 
periodic spatial grid with a Fourier pseud o-spectral solver using the s tandard 2/3 rule to avoid aliasing 
due to its quadratic nonlinearity (see, e.g.. IPatterson and Orszaa 1971 ). The pressure term is obtained by 
solving the Poisson equation 

V 2 p = -V • (u ■ Vu), (4) 

directly in Fourier space. The viscous term is treated implicitly via an exponential factor. 

To impose the no-slip boundary conditions (0 at the surface of the particle, we make use of an immersed 
boundary technique. It consists in introducing in the right-hand side of the Navier-Stokes equation (fT]) 
a penalty force fb{x,t), which acts as a Lagrange multiplier associated to the constraint denned by the 
boundary condition ([2]). The full problem ((TJ-© can then be rewritten as 

d t u = C{u) + f b , V ■ u = 0, for x e fl. (5) 

This system then has to be associated with the problem 

find fb such that fb(x,t) = for all x G r2\f2 p (t) and u(x,t) = v p (x,t) for all x 6 O p (t). (6) 

One can easily see that, when the solution is approximated by the Galerkin truncation associated with the 
Fourier spectral method, these two conditions cannot be satisfied in an exact manner: an non-zero function 
that vanishes on a full sub-domain is necessarily non-analytic and cannot be represented by a finite number 
of Fourier modes. The idea of the method which is used here is to relax the problem of finding the penalty 
force fb to 

find fb such that / \f b \ 2 d 3 x minimal and / \u — v p \ 2 d 3 x minimal. (7) 

Jn\n p (t) Jn p {t) 

All the difficulty of integrating this system is now in finding an efficient manner of solving the optimization 
problem above, which involves local condition in the physical space, while spectral methods try to be as 
local as possible in Fourier space. As this has to be done at every sub-time steps, the use of standard 
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minimization techniques would introduce a very large number of operations and result in a cost-inefficient 
algorithm. 

Efficient strategies consist in guessing acceptab le values of the field fb- We here make use of the "direct- 
forcing" method proposed by iFadlun et al. (2000). To explain this method, let us consider a simple Euler- 
discretization in time of (|5|) 

u n+1 = u n + At [C{u n ) + /"] . (8) 

The idea of the direct-forcing method is to use the fact that we want to impose u n+1 = v™ +1 in tt p to write 
from © 

f h n = -%,(*) [C{u n ) + (u n v n+1 )/At] , (9) 

where Xp( x ) denotes the characteristic function of the particle: Xp( x ) = 1 if as G f2 p (t) and otherwise. 
Other similar approaches c onsist in imposing the boundary conditions by means of a Darcy term, i.e. 



Xp( x ) ( u — v p ) (see lAngot et all I1999I ). The particle volume is then seen as a porous media where 



a large value of the parameter a ensures the no-slip boundary conditions but results in choosing a sufficiently 
small value of the time step. The advantage of direct forcing is that it consists in combining the value of 
the parameter a that is optimal for a given choice of the time step At, together with a compensation of the 
overall evolution of the velocity field inside the particle due to £. 



2.1. The smoothed particle 

The expression © for fb would be correct if the position of the grid points coincided with the immersed 
boundary. This is not true in the case of a spherical particle embedded in a cubic grid. An interpolation 
procedure is then needed in order to prescribe the bou ndary conditions at t he gr id-points next to the su rface. 
For this we use a volume- fraction scheme proposed bv lFadlun et al. ( 2000h and Pasauetti et al. ( 2008 ). The 
characteristic function Xp describing the geometry of the particle is averaged over two grid cell: 



Xp( x ) 



(2hy 



— h J —h J —h 



X P (x + y) d 3 y. 



(10) 



This representation has two advantages. First, as a typical velocity field around a particle is continuous but 
not differentiable at the surface, Gibbs oscillations are created in the vicinity of the boundary; replacing 
the characteristic function by Xp reduces such oscillations. Second, this volume averaged Xp allows for a 
continuous displacement of the immersed boundary which is important in the case of a moving particle. 



2.2. The pressure predictor 

Another issue arises when combining a spectral with a penalty method. Although applying the force 
fb to the velocity field imposes the boundary conditions precisely, the divergence-freeness of u requires 
a subsequent projection of the velocity field onto its solenoidal part which violates the formally correct 
boundary conditions. A solution to this conflict is to s olve the Poisson equation ^ for pressure with 
modified boundary conditions as for example proposed bv lTaira and Coloniusl (|2007t ) in combination with a 
finite-difference scheme. However, this is not possible with a Fourier-spectral method as in this case ((H) is 
solved with generic periodic boundary conditions. 

We attenuate this intrinsic problem of a Fourier -spectral approach by using a predictor for the pressure 
gradient in (0, as considered bv I Brown et al. ( 2001 ). We will now explain our procedure. The usual strategy 
in spectral codes without immersed boundaries consists in advancing first the velocity without taking into 
account the pressure gradient in (fTJ) 



u n+1 - 



= u n + At [-(u n ■ V)w" - vAu n ]. 



(11) 



The intermediate velocity field u n+ is compressible and is projected in a second step onto its solenoidal 
part 



n+1 



\7p n . 



(12) 
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Figure 1: (Left) Profile of the velocity (normalized to the inflow speed) in the direction of the mean flow on the axis of symmetry. 
The solid line corresponds to a computation with predictor p n ~ 1 while the points correspond to a computation without the 
predictor. (Right) L2 norm of the velocity inside the particle as a function of the number of grid points on its diameter 



This projection is equivalent to subtracting the gradient of the pressure that is given by the Poisson equation 
(TJ}. This splitting of the integration is called fractional-step method. 

If the domain contains immersed boundaries with no-slip conditions, the penalty force has to be 
additionally applied. This can be seen as an intermediate step 

u n+1 = B(u n+1 ) =u n+1 + f[ l (13) 

leading to a compressible ii n+1 , which is again projected onto its solenoidal part as in (fT2")l . A problem 
arises since the two operators P and B do not commute. The application of P spoils the result of B. Our 
approach consists in attenuating this conflict by adding a predictor for the pressure gradient to (ITT]) . We 
make of Vp™ -1 , i.e. of the value of the pressure gradient at the last time step. This predictor already 
includes the major part of the pressure needed to keep the velocity field very close to being divergence-free. 
The predictor is updated during the projection step according to 

p n = V n - X + $™, where V 2 $™ = V • ii n+1 /At. (14) 

It is important to stress that we apply the operators P and B and also the update of the pressure predictor 
during each sub-time step of the third-order temporal scheme (|3"j). 

The importance of the predictor for a Fourier-spectral scheme is evident from Fig. [TJ where the streamwise 
velocity profile on a line passing through the particle is shown. The flow conditions correspond to a numerical 
wind tunnel experiment where a fixed spherical particle (v p — 0) is embedded in a homogeneous flow. The 
velocity fluctuations introduced by the wake of the particle are eliminated in a zone at the right boundary 
by an additional application of the penalty method. The accuracy of the upstream boundary condition 
at the particle surface is strongly affected by the pressure predictor. With the predictor we achieve a no- 
slip condition at the surface while the fluid significantly enters the particle when we neglect the pressure 
predictor. Also the smoothing of the velocity fluctuation at the exit of the wind tunnel is much more efficient 
with the pressure predictor. 

A measure for the accuracy of the boundary conditions is given by the L2-norm 

v i/a 

/ u{xfdx\ (15) 

of the velocity inside the fixed spherical particle. This norm is represented in Fig. [TJ (Right) for a flow with 
a particle Reynolds number Re p = 20 and as a function of the number of grid-points ^ p along the diameter 
of the sphere. This number quantifies the resolution of the geometry of the particle by the underlying 
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homogeneous grid. The norm u 1 ^ 2 shows a power-law decrease with an exponent —3/2. This value can be 
understood by remembering that a function f(x) in a Fourier scheme is represented by a Jk{x), a truncated 
series of K modes. Using a function f(x) — max(0, cos(x)) which resembles qualitatively the flow structure 
in the vicinity of a particle at rest one can show that J B f^(x)dx cx K~ 3 . 

2.3. The forces on the sphere 

The force that is exerted by the fluid on the particle is given by the integral over the surface of the 
particle of the fluid stress tensor 

F= V-Td 3 i= / T-dS, where T = -pl 3 + v (Vu + Vu T ). (16) 

JQ p (t) JdQ p (t) 

This total force F can easily be computed from the penalty force /& because 

/ V-Td 3 i= / f b d 3 x. 

In order to distinguish the pressure and viscous terms appearing in the stress tensor (|16p. our strategy 
consists in computing the pressure contribution by constructing a homogeneous grid of points at the sphere 
surface. Then, we integrate the directed pressure pn, n being the normal vector field at the surface of the 
particle, over this grid. The viscous part is obtained by subtracting the pressure contribution from the total 
force. 



3. Benchmarks 

Now we are going to present several types of benchmarks of the proposed scheme. For this we turn 
back to the numerical wind tunnel experiment already mentioned before. The principle flow configuration 
is depicted in Fig. [2] above the critical Reynolds number where a non-stationary wake develops. 




Figure 2: Turbulent wake at Re=400 (volume rendering of high-vorticity isosurfaces). 



3.1. Drag coefficient 

An important quantity of the dynamics of a solid sphere exposed to a homogeneous inflow is the drag- 
coefficient 

where Fd denotes the total force exerted on the particle in the stream-wise direction, U c the inflow velocity, 
d is the sphere diameter, and pi the fluid mass density (set to one in our numerical settings). The measured 
coefficients are given in Fig. [3] in comparison to the empirical formula 

c ™ = rT ( 1 + - 15i?e P 687 ) ( 18 ) 

proposed by Schiller and Naumannl (jl933h which fits well experimentally measured coefficients. For all 
Reynolds numbers considered (from Re p = 2 to Re = 400) the DNS coefficients are in good agreement with 
the fitting formula. 
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Figure 3: (Left) Drag-coefficients Cd as a function of the particle Reynolds number; in parenthesis is given the number of 
internal grid points. Inset: ratio between the viscous contribution Cp to the drag coefficient and the contribution Cp from 
pressure. (Right) Convergence rate as a function of the number of grid points # p inside the sphere. 



The numbers in parenthesis in Fig. [3] (Left) denote the number of grid-points # p on the diameter of 
the sphere. The higher this number the better is the resolution of the boundary layer around the sphere. 
The thickness of this layer decreases as i?e -1 / 2 thus the higher the Reynolds number the more points =#= p 
are required in order to resolve the boundary layer. In this experiment the sphere is fixed and the inflow 
velocity is one. The Reynolds number is varied by changing the viscosity of the fluid. 

In order to check to performance of the scheme in case where the particle is moving we changed the 
frame of reference of this experiment to the particle frame and thus move the particle in a flow at rest. 
Also the zone where we remove the velocity fluctuations is moved according to the particle position. The 
measured drag-coefficients and the velocity profile around the sphere are nearly indistinguishable from the 
former wind tunnel experiments. We observe fluctuating deviations from the values given in Fig. [3] of less 
than one percent for Reynolds number up to 50 and less than five percent up to 400. 

The inset of Fig. [3] (Left) shows the ratio of the viscous Cf to the pressure Cp contri bution to the complet e 
drag Cd — Cf + Cp . The measured ratios are in reasonable agreement with results of Le-Clair et al.l ( 1970l) . 
Here we stress that the pressure contribution is not needed for the integration of the particle dynamics and 
might only be of interest for post-processing purposes. 

The convergence of the drag-coefficient can be quantified by looking at the relative error (Cd — C SN ) /C SN 
as a function of the internal points # p . From Fig. [3] (Right) we observe that the coefficients converge with an 
exponent of -3/2 for the two Reynolds numbers considered. This is the same rate of convergence as for the 
L2 norm u B 2 of the velocity inside the particle (see previous section) which suggests that the convergence- 
rate of the proposed scheme is determined by the convergence-rate of the truncated representation of the 
C°-velocity field. 



3.2. Lift- coefficient 

If the inflow velocity of the wind tunnel is not homogeneous but possesses a constant gradient 

U y = U c + g (y - y ) 



(19) 



this linear shear flow creates a lift force Fl = F ■ e y on the particle perpendicular to the mean flow. The 
rate of shear is usually given in terms of a non-dimensional shear rate 



U c 



(20) 



where g is the dimensional shear in (fT9| . Recently, the corresponding lift-coefficients 



(21) 



where Fl represents the amplitude of the lift force, have been measured by thre e different groups (see 
Dandy and Dwver , 199fi Kurose and Komori , 19991 Bagchi and Balachandarl [2002 ). 
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Figure 4: (Left) Lift-coefficients Cl as a function of the particle Reynolds number for s = 0.2. (Right) Stationary rotation rate 
(normalized to the ambient rotation rate Qf = g/2) as a function of the particle Reynolds number denned with the mean fluid 
velocity gradient ReQ = gd 2 jv. Inset: Logarithmic representation. For the exponential fit we used a = 0.8 and b = 0.043. 

From Fig [4] (Left) one observes a quite large spread of the lift coefficients measured by the three groups. 
One reason for these differing values is the sensitivity of the lif t force on the stresses in the flow and thus on 
the size of the com putational domain. It has been claimed bv lBagchi and Balachancfar ( 2002 ) that the box 
size in the work of Dandy and Dwver |l990) was too small, leading to an overestimation of the lift force. 
We have computed the lift force for two different box sizes and observe a decreasing force with increasing 
box-size. Our measurements are in good agreement with the two recent works and reproduce in particular 
the transition to negative values of the lift force for Re p « 60 — 80. 



3.3. Torque- free rotation rate 

As a last benchmark we present in Fig. @] (Right) the torque-free rotation rate in a linear shear flow 
described in the precedent section. The sphere is fixed but allowed to rotate freely. After a relaxation 
period we observe a rotation rate which is exponentially (see i nset of Fig. 0] ( Right)) decreasi ng with the 
Reynolds number. The measured values agree well with those of Bagchi and Balachandarl ( 2002 ). The small 
deviations might again be a result of different box-sizes used. 



4. Statistical properties of the turbulent wake at moderate Reynolds numbers 

This section is dedicated to the study of turbulent velocity statistics in the wake of a spherical obstacle. 
The development of a Fourier-spectral code that integrates the fluid motion and prescribes no-slip boundary 
conditions on the body surface allows one to achieve precise direct numerical simulations of turbulent wakes 
at a reasonable computational time. 

4-1. Mean flow and wake similarity 

Most models for turbu l ent fl ows past an obstacle rely on studying solutions to the Reynolds-averaged 
equations (see, e.g., IPopeL 120001) . This approach consists in writing the velocity field u(x,t) as a sum of 
its stationary average U(x) = (u(x,t)) and turbulent fluctuations u'(x,t). The mean velocity is then a 
solution to 

(U ■ V) U = -V • t - V(p) + isV 2 U, (22) 
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where r designates the Reynolds stress tensor Tij(x) = (ti£(sc,4) u'-(a;,t)). In order to obtain the properties 



of the mean velocity profile U, Eq. (|22|l is generally closed by means of the eddy viscosity assumption. This 
hypotheses consists in relating linearly the anisotropic part of the Reynolds stress to the mean rate, that is 



vt{x) [diUj + djUi 



(23) 



where k = (1/2) tracer = (1/2) (u^ u^) is the turbulent kinetic energy and the positive quantity vr{x) is the 
eddy viscosity. When the eddy viscosity z^x is assumed to be constant in space, and in particular independent 
of the streamwise distance x from the particle center, the Reynolds-averaged equation (|22p can be integrated 
to o btain the fo l lowin g self-similar estimate for the streamwise mean velocity U(x,y,z) = (u x (x,t)) (see, 
e.g- ISchuchtingL Il979h 



U(x,y,z) = U, 
where the wake half- width 7*1/2 is given by 



C D d 2 



32 vt(x — xq) 



exp 



y 2 + z 2 
r 2 ^ 2 (x) ^ 



ri/i{x) = d y/2U c (x - x )/(it In 2). 



(24) 



(25) 



Such formula, which are very commonly used in industrial applications, have two important consequences. 
First, the mean velocity deficit AU = 1 — U(x,0,0)/U c should behave at large distances in the turbulent 
wake as a; -1 . Secondly, the normalized velocity deficit profile f(x, y, z) = [U c — U(x, y, z)]/[U c — U(x, 0, 0)] 
should collapse to the similar profile f(x,y, z) = 2~^ where £ = r/ri/ 2 (x) and r — \/y 2 + z 2 . 

An analysis of the turbulent wake from a direct simulation performed at the moderate value of the 
particle Reynolds number Re p = 400 shows that this similar eddy-viscosity theory is able to rightly predict 
the decay of the mean velocity deficit. Figure[S](Left) represents the behavior of AU as a function oix/d—X. 
Its displays a power-law decay with a exponent —1 over nearly half a decade. The value —1 is undoubtedly 
obtained as seen from the inset of Fig. [5j which shows the value of the local slope (logarithmic derivative) 
of the velocity deficit as a function of the distance in the wake. The oscillations around the value —1 are a 
signature of the Gibbs phenomenon induced by the strip forcing at x w 12d 
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Figure 5: (Left) mean velocity deficit AU = 1 — U(x, 0, Q)/U c in the axis of symmetry y = and z = of the wake as 
a function of the normalized streamwise distance x/d. As seen from the inset, which represents its logarithmic derivative 
[din AU]/[dln(x/d — 1)], it displays a power-law behavior ocl/(x — xo) with xo ~ d. (Right) mean velocity deficit profile / as 
a function of the reduced spanwise distance £ = r/r\/2 for various values of x/d as labeled; the bold dashed line represents the 
Gaussian functional form of / obtained from a constant-eddy-viscosity similar profile. 



However, as seen from the right-hand panel of Fig. [SJ the mean velocity deficit profile does not seem 
to display any similarity. The half width r 1 / 2 was estimated here as the value of the distance to the 
axisymmetric axis of the wake such that the velocity deficit is exactly equal to 1/2. By definition, the 
various curves should go through / = 1 at £ = and / = l/2at£ = l. Despite these two constraints, 
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it is seen that they do not collapse on the top of each other for values of £ close to 1/2 and 3/2. The 
discrepancy is even stronger for £ > 2, where / tends to negative values that depend on x. These deviation 
at large £'s might be due to the periodicity of the spatial domain in the y and z directions that results 
in a mean streamwise velocity of the fluid that is larger than U c at y = ±7r and z = ±7T. However, the 
de viation from the simila r profi le at moderate value s of £ a re comparable to those obtained experimentally 
by|Ub eroi and Frevmuthl (|l970f ) and IWu and Faetbl (|l994j ). In particular, a striking observation is that in 



all cases with moderately large Reynolds number, the mean velocity defect profile seems to vanish for £ « 2. 
Such a feature is clearly reproduced in our simulation, as seen from Fig. [5] (Right) where the various curves 
are pinched together before reaching negative values. 

4-2. Two-point statistics and turbulence decay 

We have seen above that the self-similar wake approach reproduces well the decay oc x^ 1 of the mean 
velocity deficit. However, this approach cannot be fully true as no similarity is observed in the wake. In 
order to get a better handling, we relate in this section this observed behavior to the decay of turbulent 
kinetic energy in the wake of the spherical particle. 



It is known since the seminal work of iKolmogorovl (|1941r) that the law of decay of turbulent kinetic 



energy is dictated by the princi ple ot p ermanence ot large eddies . We tollow here the pnenomenological 
presentation of this law made bv lFrischn 1995t ). If we assume that the initial turbulent velocity correlations 
follow a infrared asymptotic self-similarity of the form (u(x + r)u(x)) ~ r at large scales r, the decay 
behavior of all turbulent quantities is dictated by the value of the exponent h. Indeed, by continuity of 
velocity statistics, the integral-scale velocity should then itself behave as U — fc 1 / 2 ~ L~ h , where k is the 
turbulent kinetic energy and L is the integral scale. However, we know that the energy decays with a rate 
that can be related to velocity by dk/dt = — e ~ —k 3 / 2 /L. Replacing k by L~ 2h in this equation, we 
obtain that dL/dt ~ L~ h , so that L ~ t x ^ 1+h K This evolution of the integral scale implies that the kinetic 
turbulent energy behave as fc ~ £-2/i/( 1 +' 1 ) 

When considering the decay of wake turbulence, all the turbulent quantities have to be of course evaluated 
in planes perpendicular to the axisymmetric axis and time has to be replaced by x/U c by invoking Taylor 
hypothesis. The above considerations hence imply that if the velocity correlations (u(r) u(0)) behave as r 
in a plane x — xq located in the developing wake, then at any x > xq, the large-scale correlations should 
remember this initial form and the turbulent kinetic energy is expected to decay 

2h/(i+h)_ Figure[5] 

represents the correlation of the radial component of the fluctuating velocity at various values of x. As seen 




Figure 6: Correlation of the velocity radial fluctuation as a function of the distance r to the axisymmetric axis for various 
values of the streamwise distance x as labeled. Inset: same in semi-logarithmic coordinates showing that the various correlation 
fall on the top of each other and follow an exponential decay. 

from this figure, velocity correlations collapse at large scales on the top of each other onto an exponential 
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Figure 7: Profile of the mean turbulent kinetic energy k in the wake of a spherical particle with Reynolds number Re p = 400. 
The top figure represents the turbulent kinetic energy k in the (x, r) plane (high values are shown in red and small values in 
blue). The bottom figure shows the behavior of k along the r = axis, clearly displaying a behavior oc x~ 2 

behavior oc exp(— r/ro) that is independent on x. Such a fast decay of the correlations is associated in 
the above phenomenological formalism to a value of the exponent h = oo (an exponential is faster than 
any power law), leading to a decay of the turbulent kinetic energy k ~ x~ 2 . This behavior is confirmed 
numerically, as can be seen from Fig. [7J 

Finally, to relate such a behavior of the turbulence decay to the way the velocity deficit decreases as a 
function of the distance downstream of the particle, one has to turn back to the Reynolds-averaged equation 
(22). The mean velocity deficit AU(x) = 1 — U(x,0,0)/U c approximatively obeys the equation 

U 2 d x AU ~ d x t xx + d y T xy + 8 z t xz ~ k/L oc x~ 2 , so that AU ~ x^ 1 . (26) 

r denotes here the Reynolds tensor introduced in previous subsections. This shows that the scaling of the 
velocity deficit is consistent with that of the decay of turbulent kinetic energy. 

5. Conclusion 

The work presented in this article aims at validating a new method for integrating the incompressible 
Navier-Stokes equation with no-slip boundary conditions on the surface of a moving particle. Simulations 
that were performed with this method in the case of a spherical particle embedded in a homogeneous flow, 
show that such a strategy reproduces with a high accuracy known results on the variation of the drag and lift 
forces as a function of the particle Reynolds number and give a precise handling on the turbulent fluctuations 
that arise in the wake downstream the sphere. 

Such benchmarking will now allow us confidently attacking the physical problems that motivated the 
development of such a code. The first question that is subject to ongoing work concerns the dynamics 
of neutrally buoyant finite-size particles that are transported by a developed turbulent flow in the case 
when the diameter of such particles falls inside the turbulent inertial range of lengthscales. The drag, lift, 
slip velocity associated to such settings are still not well defined because of the fluctuating nature of the 
fluid surrounding environment. This situation still requires an important modeling effort and a first step 
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consists in understanding how the fluid velocity statistical properties are modified by the presence of the 
particle. Preliminary results confirm that a turbulent boundary layer develops around the particle and that 
its thickness is comparable to the particle size. 

An important open question that is worth mentioning here relates to possible improvements of the order 
of the method that is presented here. As we have seen, the error of the proposed method is of the order 
of Ax 3 / 2 . To go beyond this first order, a natural idea is to make use of the variational formulation that 
is described in Sec. [5] and which consists in expressing the penalty force as a solution to a minimization 
problem. To avoid Gibbs oscillations that, as we showed, are responsible for the observed error, it could be 
envisaged to relax such minimality to the surface instead of the full particle volume, allowing for instance 
the fluid flow to develop inside the particle, without affecting the no-slip boundary condition. Such ideas 
lead to new difficulties related to the non-locality of pressure and which still need to be cautiously worked 
out. 

Acknowledgments 

This work benefited from fruitful discussions with J. Dreher. Part of this research was supported by the 
French Agence Nationale de la Recherche under grant No. BLAN07-1_192604. The research leading to these 
results has received funding from the European Research Council under the European Community's Seventh 
Framework Program (FP7/2007-2013, Grant Agreement no. 240579). Access to the IBM BlueGene/P 
computer JUGENE at the FZ Julich was made available through project HB022. Part of the computations 
were performed on the mesocentre de calcul SIGAMM and using HPC resources from GENCI-IDRIS (Grant 
2009-i2010026174). 

References 

Angot, P., Bruneau, C, Fabrie, P., 1999. A penalization method to take into account obstacles in incompressible viscous flows. 

Numerische Mathematik 81, 497-520. 
Bagchi, P., Balachandar, S., 2002. Effect of free rotation on the motion of a solid sphere in linear shear flow at moderate Re. 

Phys. Fluids 14, 2719-2737. 

Brown, D. L., Cortez, R., Minion, M. L., 2001. Accurate Projection Methods for the Incompressible Navier-Stokes Equations. 

Journal of Computational Physics 168, 464-499. 
Burton, T., Eaton, J. K., 2005. Fully resolved simulations of particle-turbulence interaction. Journal of Fluid Mechanics 545, 

67-111. 

Calzavarini, E., Volk, R., Bourgoin, M., Leveque, E., Pinton, J., F.Toschi, 2009. Acceleration statistics of finite-sized particles 
in turbulent flow: the role of faxen forces. J. Fluid Mech. 630, 179-189. 

Dandy, D., Dwyer, H., 1990. A sphere in shear flow at finite Reynolds number: effect of shear on particle lift, drag, and heat 
transfer. J. Fluid Mech. 216, 381-410. 

de Pater, I., Lissauer, J., 2001. Planetary sciences. Cambridge University Press, Cambridge, United Kingdom. 

Fadlun, E. A., Verzicco, R., Orlandi, P., Mohd, J., 2000. Combined Immersed-Boundary Finite-Difference Methods for Three- 
Dimensional Complex Flow Simulations. Journal of Computational Physics 60, 35-60. 

Frisch, U., 1995. Turbulence: The legacy of A.N. Kolmogorov. Cambridge University Press, Cambridge, United Kingdom. 

Homann, H., Bee, J., 2010. Finite-size effects in the dynamics of neutrally buoyant particles in turbulent flow. J. Fluid Mech. 
651, 81-91. 

Kolmogorov, A., 1941. On degeneration of isotropic turbulence in an incompressible viscous liquid. Dokl. Akad. Nauk SSSR 
31, 538-540. 

Kurose, R,., Komori, S., 1999. Drag and lift forces on a rotating sphere in a linear shear flow. Journal of Fluid Mechanics 384, 
183-206. 

Le-Clair, B. P., Hamielec, A. E., Pruppacher, H., 1970. A Numerical Study of the Drag on a Sphere at Low and Intermediate 

Reynolds Numbers. Journal of Atmospheric 27, 308—315. 
Naso, A., Prosperetti, A., 2010. The interaction between a solid particle and a turbulent flow. New Journal of Physics 12, 

033040. 

Pasquetti, R,., Bwemba, R., Cousin, L., 2008. A pseudo-penalization method for high Reynolds number unsteady flows. Appl. 
Numer. Math. 58, 946-954. 

Patterson, C, Orszag, S., 1971. Spectral calculation of isotropic turbulence: efficient removal of aliasing interaction. Phys. 
Fluids 14, 2538-2541. 

Pope, S., 2000. Turbulent Flows. Cambridge University Press, Cambridge, United Kingdom. 

Prosperetti, A., Oguz, H., 2001. Physalis: A new o(n) method for the numerical simulation of disperse systems, part i: Potential 

flow of spheres. J. Comput. Phys. 167, 196-216. 
Qureshi, N. M., Bourgoin, M., Baudet, C, Cartellier, A., Gagne, Y., 2007. Turbulent transport of material particles: an 

experimental study of finite size effects. Physical Review Letters 99, 184502—184506. 



12 



Schiller, L., Naumann, A., 1933. fiber die grundlegenden berechnungen bei der schwcrkraftaufbereitung. Verein Deutschcr 

Ingcnicurc 77, 318-320. 
Schlichting, H., 1979. Boundary-Layer Theory. McGraw-Hill, New York, United States. 

Seinfeld, J., Pandis, S., 1998. Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. Wiley-Interscience, 
New York, United States. 

Shaw, R., 2003. Particle-turbulence interactions in atmospheric clouds. Ann. Rev. Fluid Mcch. 35, 128-277. 
Shu, C, Osher, S., 1988. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J. Comput. Phys. 77, 
439-471. 

Taira, K., Colonius, T., 2007. The immersed boundary method: A projection approach. Journal of Computational Physics 225, 
2118-2137. 

Uberoi, M., Freymuth, P., 1970. Turbulence energy balance and spectra of the axisymmetric wake. Phys. Fluid 13, 2205-2210. 
Volk, R., Calzavarini, E., Leveque, E., Pinton, J.-F., 2010. Dynamics of incrtial particles in a turbulent von Karman flow. 

accepted for publication in J. Fluid Mech., 10. 
Wu, J.-S., Faeth, C, 1994. Sphere wakes at moderate Reynolds numbers in a turbulent environment. AIAA Journal 32, 535—541. 
Xu, H., Bodcnschatz, E., 2008. Motion of incrtial particles with size larger than Kolmogorov scale in turbulent flows. Physica 

D: Nonlinear Phenomena 237, 2095-2100. 



13 



